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Abstract 



The abundances of elements in the Earth and the terrestrial planets provide the initial conditions for life 
and clues as to the history and formation of the Solar System. We follow the pioneering work of |Bond et al.| 



(2010a) and combine circumstellar disk models, chemical equilibrium calculations and dynamical simulations 
of planet formation to study the bulk composition of rocky planets. We use condensation sequence calcula- 
tions to estimate the initial abundance of solids in the circumstellar disk with properties determined from 
time dependent theoretical models. We combine this with dynamical simulations of planetesimal growth 
that trace the solids during the planet formation process and which include the effects of gravitational and 
hydrodynamical mixing. We calculate the elemental abundances in the resulting rocky planets and explore 
how these vary with the choice of disk model and the initial conditions within the Solar Nebula. 

Although certain characteristics of the terrestrial planets in the Solar System could be reproduced, none of 
our models could reproduce the abundance properties of all the planets. We found that the choice of the 
initial planetesimal disk mass and of the disk model has a significant effect on composition gradients. Disk 
models that give higher pressure and temperature result in larger variations in the bulk chemical compositions 
of the resulting planets due to inhomogeneities in the element abundance profiles and due to the different 
source regions of the planets in the dynamical simulations. We observed a trend that massive planets and 
planets with relatively small semi-major axes are more sensitive to these variations than smaller planets 
at larger radial distance. Only these large variations in the simulated chemical abundances can potentially 
explain the diverse bulk composition of the Solar System planets, whereas Mercury's bulk composition can 
not be reproduced in our approach. 

Keywords: Abundances, interiors, Planetary formation, Terrestrial planets 



1. Introduction 



The process of planet formation is far from understood. The established scenario for the formation of the 
Earth and other terrestrial planets is that most of their masses were built up through the gravitational 
collisions and interactions of smaller bodies (Chamberlin 1905 Safranov 1969 Lissauer 19931. These 
smaller bodies, often called planetesimals, form via the coagulation of dust in the midplane of a circumstellar 
disk. This transition from dust particles to gravitationally interacting bodies is an important topic of ongoing 
research (Weidenschilling 1980 Johansen et al. 2007 Chiang and Youdin 2010). 

Many direct simulations of the formation of rocky planets from planetesimals and protoplanets have been 



performed ( Chambers and Wetherill 1998 Agnor et al. 1999 Raymond et al. 2004 O'Brien et al. 2006 
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Kokubo et al. 2006 1 . One of the goals of these numerical studies is to reproduce the inner part of the Solar 
System, primarily the number of planets and their mass distribution and orbital parameters. Additional 



constraints can be the delivery of water rich material onto the planets ( Raymond et al. , 2004 ) or the timing 



of a Moon forming impact taking place ( Touboul et al. 2007 ) . Dynamical simulations by Morishima et al. 



(2010) include gas drag and type-I migration and the secular perturbative effects of Jupiter and Saturn. 



Some of the simulated planetary systems partially satisfy the main constraints given by our Solar System 
but there are still some issues to resolve. For example, the small mass of Mars is not reproduced in most 
of the dynamical simulations, although a possible solution to this was recently proposed by |Walsh et al.| 



Determining the composition of the terrestrial planets is a challenging task. In case of the Earth, samples 
from the crust can be easily obtained and lava provides samples from the mantle. There are samples from 
Mars in the form of meteorites and in situ analyses of rocks from spacecrafts and probes on Venus and 
Mars. Remote sensing and geophysical measurements provide additional data for all planets. Modelling has 
to include many different processes that are involved in creating the final bulk chemical composition of a 
terrestrial planet. It reveals that the bulk composition of the planets shows some significant differences in 
the relative abundance of the major rock forming elements (Morgan and Anders| |1980 |Kargel and Lewis 



1993 Lodders and Fegley 1997) or their metallic core masses, e.g. Righter et al. (2006). Most meteorites 



and the rocky planets show a depletion in volatile elements relative to primitive carbonaceous CI meteorites, 
which resulted from the formation of solids in the protoplanetary disk ( |Davis| |2006[ ) or from volatile loss 
during the planet formation process. Albarede et al. (2009) pointed out that the volatile depletion of the 



Earth can be explained by a very hot protoplanetary disk in the inner part of the Solar System, which 
suppresses the condensation of volatile elements. On the other hand, the Earth is enriched in refractory 



elements (Palme and O'Neill 2003 Rubie et al. 2011). 



When considering the final elemental compositions of the planets it is also important to consider the effects 
of dynamical evolution. The planetesimals that eventually come together to form a terrestrial planet, 
themselves formed in various locations in the protoplanetary disk. Their bulk composition can depend 
strongly on their initial location since the formation of solid species in the protoplanetary disk will most 
likely depend on the local gas composition, temperature and pressure. Gravitational scattering among the 
planetesimals and protoplanets and radial movement due to interaction with the gas disk results in a radial 
mixing of the building blocks of the planets and in different bulk chemical compositions of the final planets. 



Bond et al. (2010a) introduced the bulk chemical composition of the Solar System's planets as an additional 



constraint on the models of their formation. Based on an a-disk model by Hersant et al. (2001), they 



determined the spatial equilibrium composition of solid material in the inner protoplanetary disk with the 
HSC Chemistry software package. These spatial solid material abundances were mapped onto the initial 



planetesimals in the dynamical simulations of O'Brien et al. (2006). The initial location of each planetesimal 



that is embodied in the final planets is traced in those simulations. Assuming that the bulk composition of 
the final proto-planets is the sum of the solids embodied in the initial planetesimals, they estimated the bulk 
chemical composition of the simulated planets. They found that these are comparable to the bulk chemical 
composition in the Solar System in the case of most elements, although volatiles are an exception and were 
often overestimated. Furthermore, Bond et al. (2010b) applied their method to extrasolar giant planet 
systems. They performed dynamical simulations designed to study the formation of rocky planets around 
those stars with known gas giants detected via radial velocity and known abundances from high resolution 
stellar spectra. They estimated the bulk chemical composition of the simulated planets and sometimes 
extreme variations in composition were found. 



Bond et al. (2010a 



In this report, we followed the methodology of Bond et al. (2010a). We combined different disk models and 



the dynamical simulations performed by Morishima et al. (2010p to estimate the bulk chemical composition 
of the simulated planets. 



were interested in "making the Earth" and reproducing 
abundances trends, geochemical ratios and oxidation states of the planets. In addition, they studied the late 
veneer and the loss of volatiles in collisions. In comparison, we focused on the bulk chemical compositions 
and on variations that might result from different initial conditions in the dynamical simulations and our 
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assumptions in the modeling of the gas disk. In addition, we did not treat the transition from the disk 
model to the dynamical simulations as a free parameter. Instead, we chose a more self-consistent transition. 
With those differences, we hope to improve the methodology. In the following, our disk models give a time 
dependent description of the density, temperature and pressure profile of a protoplanetary disk - the physical 
conditions for the equilibrium chemistry calculations. 

The goal of this work consisted of three parts: First, we wished to infer how the distribution of solids in the 



radial direction in the protoplanetary disk depended on the choice of the disk model. Bond et al. (2010a I is 
only based on one model, whereas we explored three simple but widely used disk models and their resulting 
chemical profiles. This allowed us to estimate how much the bulk composition of planets depends on the 
conditions in the disk. 

Second, we focused on the dynamical simulations and tried to understand the link between the initial 
locations of the planetesimals embodied in simulated planets and the planet's final position and mass. A 
detailed understanding of this link gives deeper insights into planet formation from the dynamical point of 
view. We compared the cumulated initial distribution of groups of planets to infer if this can be another 
indicator of trends in the composition of planets. 

Finally, we studied the dependence of the estimated bulk composition of the simulated planets on the disk 
model and the initial conditions of the dynamical simulations. We were especially interested in the most 
extreme cases that can be obtained, since such an extreme choice can probably have a strong effect on the 



results in bulk composition studies. Hence, we could cover a broader range of initial assumptions than Bond 



et al. (2010a). We compared our results with Solar System bulk compositions to see if this modeling can 



reproduce earlier predictions. 

There are uncertainties in the modelling of the protoplanetary disk and assumptions are made in the de- 
termination of the equilibrium composition. The transition from the disk model to the disk of solids used 
in the dynamical simulations is performed in a simple way. The dynamical simulations represent different 
sets of initial conditions and their effect on the formation process is only poorly understood. The biggest 
challenge in this work was to deal with those issues and to investigate how the final results were affected by 
them. We also keep in mind that fundamental processes like disequilibrium condensation and core-mantle 
differentiation are not taken into account although they might play a crucial role in determining the bulk 
chemical composition of a planet. 

This report is structured as follows: In section[2j we review the disk models that are used. We briefly present 
the tool to calculate the equilibrium abundance in the disk and the dynamical simulations. In section [3j the 
bulk chemical compositions of the planets are shown and the differences that arise from the different disk 
models and initial conditions. In section [4j the caveats and implications of our approach are summarized. 
Finally, the most important conclusions are given in section [5j 



2. Methodology 

Our method consists of a two-step approach first treating the condensation of solids in an evolving circum- 
stellar disk and then following the reshuffling of those materials dynamically through the planet formation 
process. The evolution of the surface density, pressure and temperature profiles in the disk are given by 
a disk model. Based on these profiles, a model that describes the chemical processes in the disk gives the 
abundance of solid species in the disk midplane in space and time. The disk can evolve until it satisfies a 
selected criterion based on the initial surface density in the dynamical simulations. At this transition time, 
the distribution of solids in the disk midplane based on the disk model gives the bulk composition of the 
initial planetesimals in the dynamical simulations. The following dynamical interactions of the particles and 
their mergers link the bulk composition of the initial planetesimals to the bulk chemical composition of the 
final planets. 
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2.1. Disk models 



Observations reveal that protoplanetary disks are not static structures. The observed lifetimes are several 
millions of years over which period mass and angular momentum are transported in the radial direction. 
The evolution of the disk's surface density £ follows from the conservation of angular momentum and mass 
and can be expressed in the evolution equation, which is given in the case of a thin disk by |Lyn dcn-Be lTand| 
Pringle | ( [T974| : 

as 3 a / 1/2 o 



Of 



3_9 
r dr 



Or 



(1) 



The angular momentum transport is controlled by the viscosity v. In the widely-used a-prescription, the 
viscosity is expressed in the form 

v = ac s H = ac^Q- 1 (2) 

where c s is the local sound speed and H is the local scale height of the disk. They are related by c s = HH, 
with the Keplerian angular motion 

Igm q 



n 



(3) 



where G is the gravitational constant and Mq is the mass of the central star, which in our model we 
always assume to be a solar mass star, a is a dimensionless quantity which controls the efficiency of 



angular momentum transport due to viscosity and was introduced by Shakura and Sunyaev (1973). This 



parametrization reflects the uncertainties about the processes causing viscosity. The energy dissipation due 
to viscosity is one of the main heat sources in the disk, especially close to the central star. At larger radii, 
the irradiation by the central star becomes dominant. 

The isothermal sound speed is 



' hT c 



(4) 



where k\, is the Boltzmann constant, /i the mean molecular weight and run the mass of a hydrogen atom. 
We assume that all the energy is dissipated close to the midplane of the disk with T c being the midplanc 
temperature. 

Two and even three-dimensional radiative transfer models of circumstellar disks are powerful tools to explore 
their physical structure and chemical evolution. Very sophisticated models of circumstellar accretion disks 
were developed in the last years (see Dullemond et al. (2007) for a review). Here, we restrict ourselves to 
simple analytic models. We explore three different disk models here: Two analytic one-dimensional models 
and a two-dimensional model. In what follows, we derive the surface density, the temperature and the 
pressure profiles that result from each model. 



2.1.1. Self similar solution 



The self-similar solution of the evolution equation ([I]) derived by Lynden-Bell and Pringle ( 1974 1 has the 
form 

S fo *) = ~vH 5 /2-7)/(2-7) exp (J^\ , (5) 



3irv 



A 



where 



A = 



1. 



ts = 



3(2-7) 2 v 



(6) 
(7) 
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with f — r/r a , where r s is the scale length. It gives the transition radius at which the disk's exponential 
cut-off starts. This self-similar solution holds if the viscosity is given by a power law, v (r) oc (r/r s ) 7 . We 
adopt 7 = 1. In the a-prescription, |2| and Q result in a temperature profile with a fixed power law in r: 

T c (r) oc f- 1 ' 2 . (8) 

Then, v(r s ) is obtained by d2j). C is a constant and can be estimated if the initial disk mass is known. 

In principle, the above equation describes a time-independent midplane temperature. Since the surface 
density evolves in time, we are interested in a time-dependent temperature profile. Hence, we use the above 
power-law of T c cx r^ 1 / 2 and derive a surface density-dependent normalization of the profile. To estimate 



the temperature at r s , we follow Armitage (20101, assuming that the disk surface emits like a black body 
and the vertical energy flux equals the dissipation rate at the disk midplane and is constant in height. The 
temperature of the disk at r s is given by 

T^r s ,t)= 2 ^^ K ,(r s Mr s ,t) 2 r^, (9) 

where k is the opacity and a is the Stefan-Boltzmann constant. The normalization of the temperature profile 
becomes time dependent. This approach is not fully self-consistent but it provides a simple model. Finally, 
the pressure P is given by the ideal gas equation of state by: 

|/5tom (10) 

2.1.2. Chambers model 

A problem that occurs in the self-similar solution model is that the viscosity is constant in time even when 
the disk cools and the opacity is constant with radius even very close to the star, where all grains will 



sublimate. Chambers (2009) describes another analytic model which takes into account heating by viscous 
accretion and radiation from the central star. Three different regimes are introduced: in the innermost disk, 
where dust grains are sublimated, the opacity follows a power law. Then, for the intermediate part of the 
disk, a constant opacity is assumed. Viscous heating is the dominant heat source in these regimes. In the 
third regime, the outer disk, stellar irradiation is the most dominant heat source. 

In order to obtain a first approximation, in most parts of the disk we use a constant opacity kq. Although 



more complex, piecewise continuous opacity laws are used widely (Ruden and Pollack (1991), Stepinski 



(1998)), there are still uncertainties in the details. 



In the intermediate regime, where only viscous heating takes place and the opacity is constant, the disk 
surface density and temperature are given as a function of time and radius as follows: 

/ \ -3/5 / , \ -57/80 

r \ -9/10 / t \ -19/40 



T(r,t)=T vis [ g -J ( 1 + ^J (12) 



S v is and T V i S are the initial surface density and temperature at the outer disk radius. They are given by 

7Af disk _(0) 
s 2 



Evis " 10ns 2 ' (13) 



27 KOjkB\ 1/3 2/3 l/3 



\64 afirriH 
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where 7 ~ 1.4 is the adiabatic index, Do = y/ GMq / Sq is the Keplerian velocity at the initial disk radius 
and so is the initial outer radius of the disk, which is the scale length of this model. 

The viscous time scale r vis is given by 

1 /J-rriH QqMq 



167T ajk B £ vis T v 



(15) 



As the surface density of a disk decreases with time, the stellar irradiation will become a dominant heating 
source and an initially fully viscously heated disk will start to become mostly heated by irradiation in the 
outer regions, beyond a radius r t . Then, at r > r t , 

\ -15/14 / . \ -19/16 

r \ I t • 

21(r,t) = L rad I 

where 



=(r, t )-H-l-l 1 + (16) 



with 



Moreover, 



S ra d — S vis ( — ^ J , (17) 
\ J- rad / 

\ -3/7 

T(r.l) I [ ) . (19) 



The radius r t can be expressed as 



•y \ 70/33 / , x -133/132 

rt(*) = -o(i^J 1 + . (20) 



r. 



VIS 



and this transition radius r t moves inwards with time. These equations hold if the initial outer edge of the 
disk is always in the viscous regime, sq < r tl which means that irradiation can be initially neglected at the 
outer disk edge. This is a suitable assumption in our case. 



In the regime very close to the star, Chambers (2009) assume that the opacity k depends on the temperature 

K = K o(^) , (21) 



Te, 

where T e is a transition temperature which defines the border between the inner most region and the 
intermediate regime. Stepinski (1998) recommended T e = 1380 K and n = —14. The surface density and 
temperature change with time as follows: 

\ -24/19 / , \ -17/16 

r \ / t x 



EM-JW^-j + , (22) 

with 

/ T n 14/19 

^cvap — ^vis I 1 7 (23) 

and 

/ \ -9/38 / , \ -1/8 

The transition radius between the inner and the outer viscous regimes has the form 

/ v \ 95/63 / , v -19/36 

re(t) = So (^£) 1 + — (25) 

\ •'-'vis / \ Tvis / 

and moves inwards with time. 



2.1.3. Two-dimensional model 



The analytic models are given by radius-dependent quantities and basically all quantities in the vertical 



direction are represented by their central or averaged values. Following Papaloizou et al. (1999), Hersant 



et al. (2001) and Hure (20001, we estimate the pressure and temperature profiles of the disk in a more 



sophisticated way in a two-dimensional model: in radial direction r and in vertical direction z. For each 
radius, a set of differential equations describes the vertical structure of the disk, where we ignore self-gravity. 



The structure is described (Armitage 2010) by hydrostatic equilibrium 



by the vertical variation of the flux F(z) 



1 dP 

p dz 



dF 9 



= -n 2 z, 



p{z)v(z)Vt 2 



dz 4 

and by the radiative heat flux in the diffusion approximation for radiative transfer, 

dT 3np(z) 



dz 16aT(z) : 



;F(z). 



In addition, there is an equation of state 

P(z) 



p{z)kBT 4cr 



T(z) 4 



(26) 



(27) 



(28) 



(29) 



pmH 3c 

The second term is the radiation pressure which is usually very small and we ignore it. Turbulent pressure 



can be ignored if a <§C 1, Hure (2000) 



The opacity n is in general a function of P and T. In this model we use an opacity model presented by 



Ruden and Pollack (1991), which is a piecewise, continuous function of temperature. More sophisticated 



models add an additional atmosphere-like optically thin layer above the optically thick disk (Hure 2000) 
which we neglect. 



Papaloizou et al. ( 1999 ) used the above set of equations to estimate the vertical structure of protoplanetary 



disks. The pressure P{H) is, with a boundary optical depth of r(if) = 2/3, given by 

2Hil 2 



P(H) = 



3 k 



(30) 



Papaloizou et al. (1999) pointed out that the results do not depend on the exact choice of t(H) if t(H) <C 1. 



At least the midplane pressure does not depend sensitively on this value, although this need not be the case 
for values at other scale-heights. 



The accretion rate M acc (t) in the disk decreases in time according to (Drouart et al. (1999)), 



M acc (t) = M acc (0) 



1 + ^ 



-3/2 



where 



R D 



(31) 



(32) 



Rd is the time dependent outer radius of the disk. The initial accretion rate is assumed to be M acc (0) 
5 x lO- 6 M yr- 1 . 
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Element Abundance Element Abundance 



H 


12.00 


Si 


7.51 


He 


10.93 


P 


5.36 


C 


8.39 


S 


7.14 


N 


7.78 


Ca 


6.31 





8.66 


Ti 


4.90 


Na 


6.17 


Cr 


5.64 


Mg 


7.53 


Fe 


7.45 


Al 


6.37 


Ni 


6.23 



Table 1: Relative element abundances in the present-day solar photosphere from |Asplund et el!.| ( |2005[ l ■ The abundance is given 
as the exponent to the base 10. Hence, H has the highest abundance in the solar photosphere, its relative abundance is 10 12 . 
If this value is used as input for the HSC Chemistry module, it can be interpreted as the number of moles of H, whereas the 
normalization is arbitrary and can be neglected. 



Each annulus of the disk can be calculated independently. To solve the set of equations ( 26 - 29 ) , initially, a 
scale height H is guessed. Calculations start at the top of the disk located at height H , where the boundary 
conditions are given by the optical depth t(H) = 2/3, and 



F(H) = — M acc (t)n 2 . 

87T 



(33) 



In the midplane, the energy flux along the vertical direction must be -F(O) = due to disk symmetry. If the 
estimated flux F est (0) in the first calculation is not sufficiently small, that is \F est (0)/F(H)\ < 10~ 4 , a new 
H is chosen, estimated through a Newton- Raphson method. Given the initial disk radius Rr>, a-parameter 
and the initial accretion rate, we can thus calculate density, temperature and pressure profiles at any time 
t. Hure and Galliano (2001 ) reports of a 30% maximum deviation on the midplane temperature comparing 



one-dimensional and two-dimensional models. 



2.2. Chemistry 



To estimate the chemical composition of solid material in the protoplanetary disk midplane, equilibrium 
condensation of the gas is assumed. In a closed system of elements, condensation calculation models de- 
scribe the equilibrium distribution of the elements between coexisting phases. For decades, condensation 
calculations have been used to describe chemistry in astrophysical problems, for example the formation of 
the Earth ( Latimer 1950 ) . Equilibrium condensation provides a basic framework to predict and interpret 



bulk chemical composition of meteorites and planets ( Ebel 



2006 Davis 


2006 


Righter et al. 


2006 


). Hence 


study. Following previous studies ( Pasek et al. 



2005 Bond et al. 2010a), we use the commercial HSC Chemistry software package to perform the equilib- 
rium condensation calculations. The equilibrium chemical composition of a system of elements is found by 
iteratively minimizing it's Gibbs free energy, which is done in the GIBBS equilibrium solver. The algorithm 
is further described in White et al.| ( 1958 1. Over one hundred different gaseous and solid species are taken 



into account, see Bond et al. (2010a). The initial parameters for the HSC Chemistry equilibrium composition 
calculation module are the amount of each pure gaseous element in kmol. It is assumed that no other species 
are present initially. In addition, the solids that are present when equilibrium is reached are also specified. 
Starting from a gas of solar composition ( Asplund et al. 2005 ) , see table [I] and under the assumption of 
chemical equilibrium, we calculate which distribution of compounds is most thermodynamically stable at 
midplane pressure P(r,t) and temperature T(r,t) given by the disk models. We do the calculations in a 
radial range of r = [0.25 AU,4 AU] with a step size of Ar = 0.05 AU. 

The output is the amount of each compound in kmol. We are interested in the relative abundances of the 
elements in the condensates. Therefore, we weight the abundance of a certain species with the number of 
atoms of the element included in that certain molecule and its atom weight. The relative abundance is 
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expressed as wt.% which means weight percent(age): 



^.%oiX= s P^ lxX ^ r ]. (34) 

where ha is the total abundance of species A located at radius r and Xa gives the number of X-atoms in 
species A. mx is the atom weight of the element. The numerator is a summation over all solid species that 
contain the element X. The denominator is a summation over the weight of all atoms Y in all species A in 
the solid phase at radius r. 

In the following, we often order the involved elements according to their volatility. A good measure of 
the volatility is the temperature at a given pressure at which 50 % of the original element has condensed 



(Lodders, 2003). In the case of O (and C), this temperature is lower than any temperature reached in our 
disk models, since most O condenses in water ice. Therefore, its volatility is measured according to the 
condensation of the higher temperature condensates. 



2.3. Dynamical simulations 



Morishima et al. (2010) carried out 64 simulations which describe the collisional growth of planetesimals and 



the subsequent long-term evolution and stability of the resulting planetary systems. The simulations explore 
sensitivity to the initial conditions, including the initial mass and radial distribution of planetesimals, the 
timescale for the dissipation of the solar nebula and different orbits of Jupiter and Saturn. The simulations 
start with 2000 equal-mass particles placed between 0.5 and 4 AU embedded in a gas disk. The planetesimal 
disk has an initial mass TOpi anet esimais of 5m e or 10 m 9 . Since we use this mass very often in this report, 
we introduce m5 for the case of a dynamical simulations with fripianetesimais = 5 rn® or mlO for the case with 



wher 



e p 



10 m® The surface density £ of solids and of the initial gas disk depends on the radius, S oc ■ 
is 1 or 2. The gas disk dissipates uniformly in space and exponentially in time with a gas dissipation time 
scale T gas = 1, 2, 3 or 5 Myr. After time r ga s from the beginning of the simulation, when a significant part of 
the gas disk has disappeared, Jupiter and Saturn are introduced on their orbits. Several orbits are tested: 



CJS are Circular orbits according to the initial conditions used in the Nice model (Tsiganis et al. 2005) 



EJS place Jupiter and Saturn on their current orbits, EEJS are the same as EJS but higher eccentricity for 
Jupiter, CJSECC are the same as CJS but higher eccentricities for Jupiter and Saturn. Each simulation 
represents one unique combination of initial conditions. See tables 1 and 2 in Morishima et al. (2010) for 
further details. 



The simulations in Morishima et al. (2010) were constructed in an attempt to reproduce the global char- 
acteristics of the terrestrial planets in the Solar System. The constraints which were used to compare the 
simulation results with the observations are the spatial mass distribution, the deviation from circular and 
coplanar orbits and the timing of the Moon- forming impact. In general, 1 — 5 planets are formed in the 
terrestrial region. An observed trade-off is that if a similar radial mass concentration is achieved, the Moon- 
forming impact occurs too early. The dependence on the initial conditions becomes manifest for example in 
the fact that small gas dissipation timescales (r = 1 — 2 Myr) are needed to avoid a significant depletion of 
solid material due to migration. 

Other simulations ( |Q'Brien et al. 2006 ) often consist of initially roughly Mars-mass embryos embedded in 
a planetesimal disk with no mutual interaction between the planetesimals. The initial mass of the fully 
interacting particles in the Morishima et al. (2010) simulations is smaller than the Moon mass. Hence, 



Mercury analogs concerning its mass and semi-major axis can form in those simulations in contrast to other 
studies. Nevertheless, these calculations assume perfect accretion and no mass loss due to giant impacts - 



such an event is most probably responsible for striping off Mercury's initial crust and mantle (see Benz et 



al. ( 2007 ) for a review) 



The Minimum Mass Solar Nebula (Hayashi, 1981) has a gas density profile £ oc r~ p with p = 1.5. The 
steady-state solution of ( lj) pr edicts p = 1. Observations suggest a flatter profile with p < 1 (Kitamura et 
al. 2002 Andrews et al. 2009 ). Therefore, a gas density profile with p = 2 is extremely steep in the context 
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of circumstellar disks and we focus on the dynamical simulations with p = 1. The initial surface density of 
solids at 1 AU is So, solid = 6.1g/cm 2 in the case of m5 or Eo )SO iid = 12.2 gcm~ 2 for mlO. The gas density at 
1 AU is S , gas = 2000 gem -2 for both m5 and mlO. 



2.4- Transition from disk models to dynamical simulations 

In order to combine the disk models and the dynamical simulations, we have to choose a suitable transition 
criterion. Since our goal is to study the evolution of the solid material in the dynamical simulations, we 
define the solid surface density at 1 AU as the normalization condition. This means that the surface density 
of solids Egoiid predicted by the disk model at 1 AU at transition time t = itrans has to reproduce the initial 
solid surface density of the dynamical simulations at 1 AU, So lSO iid- For example, 

S so iid(l AU, t d isk) = £o,soiid = 6.1 g/cm 2 , (35) 

in the case of the m5 simulations. Actually, the disk models give the surface density of gas S gas and we 
assume a gas-to-dust ratio of 

^44 = io °' ( 3e ) 

S so iid(r, t) 

which is roughly the gas-to-dust ratio of the ISM. This value is smaller than the gas-to-dust ratio in |Mor-| 



ishima et al.| ( |2010[ ); 

£ ,gas_ = 329j (37) 



J 0, solid 



for m5 simulations, and half that for mlO. These ratios are rather high and can only be achieved with a 
very high solid mass loss rate combined with a low gas dissipation rate. Since we model the solids in this 
disk, the ISM ratio seems to be a better choice. 



The initial total disk mass in the dynamical simulations is 

S 

J 0. solid 



TOtot ~ TOplanctcsimals ^ = 1650 m© = 0.0055 Mq, (38) 



for m5 and the mlO simulations. Therefore, the initial mass of 0.1 M in our disk models is a justified 
assumption since the disk loses mass with time. Disk model masses at t = itrans are shown in table [2] 

A disk model is described by the initial mass, time t, the a-parameter and the disk scale radius r s . For 
all models we adopt the following initial constraints: the initial disk mass Mdi s k(i = 0) = 0.1 M© and 



a = 0.009, Hersant et al. (20011. The models use different opacity laws. By default, the constant opacity 
is kq = 3 cm 2 g _1 . We adopt a radial scale length r s — 10 AU. We are only interested in profiles that 
reproduce the surface density at 1 AU. The choice of the scale radius has no significant effect on the results: 
In the self-similar solution model, the effect of the scale radius r s on the surface density distribution can be 
compensated by a suitable choice of itrans and C . It has a minor effect on temperature and pressure. 

In the Chambers model, the initial disk radius and the time t are free parameters and they are degenerate 
in a similar way as in the self-similar solution model. In the two-dimensional model, the accretion rate is 



a function of the initial outer disk radius, time and initial accretion rate (see equation (31)). Since the 
surface density in this model is controlled by the accretion rate, the disk model is also degenerate in the 
parameters Rd and itrans- Actually, the choice of the scale radius is not completely free: In the case of 
a very large initial scale radius (r s ~ 100 AU), the initial disk mass might be too small to reproduce the 
surface density at 1 AU for any t > 0. For r s — 10 AU, there is a itrans > such that S(l AU, idisk) = So jSO iid 
for all models. The same holds for the choice of the initial accretion rate in the two-dimensional model. 
If the initial accretion rate is too small, it will never result in a suitable surface density. Table [2] lists the 
model-dependent transition times which can differ by orders of magnitude. 



Note that in Bond et al. (2010a I, the transition time is left as a free parameter. Hence, a better fit to the 



observed planetary abundance values is possible, but this approach is not as self-consistent as ours. 
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disk model disk mass itrans[10 4 yr] MdiskflO 2 m©] r; n [AU 



self-similar solution mr = 5m® 3.7 3.5 0.50 

self-similar solution mj- = lOm® 2.1 4.4 0.75 

Chambers mr = 5m0 5.0 5.0 0.30 

Chambers mj = 10 1.9 6.0 0.50 

two-dimensional rriT = 5m^ 25.0 1.2 0.35 

two-dimensional my = lOm^ 41.0 1.8 0.65 



Table 2: The adopted transition time and the resulting disk mass for the different disk models and solid disk masses. r- ln is the 
innermost semi-major axis at which solids condense at ttrans- The model parameters are a = 0.009 and r s = 10 AU and the 
initial disk mass is m ( j; s ] I (0) = 0.1 Mq. 



3. Results 

This section presents the different disk models that are normalized to the initial conditions of the dynamical 
simulations. The resulting chemical abundance profiles are shown briefly. Then, the source regions of the 
planets are studied. Combining the chemistry and the dynamics results in the bulk chemical composition 
of the planets which are studied in detail. Two elements which are very important for the emergence and 
evolution of life, H and C, are briefly highlighted at the end of the section. 



3. 1 . Initial disk profiles 

Surface density profiles and the corresponding pressure and temperature profiles at itrans for the different 
models are plotted in figure [TJ The slope in the temperature and pressure profiles of the models are very 
similar. At radial distance r > 1 AU, the different opacity models change the slope of the temperature profile. 
The slope of the dynamical simulation surface density (E oc r" 1 ) is only reproduced in the self-similar model 
with the other models yielding flatter profiles. The high solid surface density of the mlO simulations results 
in disk models that predict a high temperature and pressure profiles relative to the m5 disk. The Chambers 
model gives the highest overall temperatures, especially at larger radii. In general, the two-dimensional 
model gives the lowest temperature. These differences allow us to explore a range of different temperature 
and pressure profiles in our study. 

Comparison with disk profiles obtained in previous studies is difficult since our initial assumptions vary. 



A brief comparison with e.g. the sophisticated study by D'Alessio et al. (1998) reveal that our profiles 
represent a reasonable approximation to more complex models. 



3.2. Elemental abundance profiles 

Figure [2] shows the elemental abundances in solids as a function of radius in the case of the two-dimensional 
model and a disk of mass m5 and the Chambers disk model in the case mlO as extrema. The two-dimensional 
model is a good case study to explain some of the main features since this model gives the largest numbers 
of solids along the disk: at 1.3 AU, troilite (FeS), at 2 AU, fayalite (Fe2SiO/i) and around 3 AU iron oxide 
(Fe3C>4) all condense out of the solar nebula. Serpentine (Mg 3 Si205(OH)4) and aluminum oxide (AI2O3) 
follow at 3.1 AU. Since all models have similar temperature profile slopes at r > 1 AU but predict different 
temperature normalizations, the location of condensation shifts to larger radii in the case of colder disks. 
In the case of the Chambers model in figure [2] the condensation of FeS is the only feature beyond ~ 1 AU 
beside the occurrence of Fe, magnesium silicates and other species around 1 AU. 

We point out that the total amounts of solid material changes completely with semi-major axis. Hence, 
traces of refractory solids that condense close to the star can dramatically change the percentage weights of 
the different elements. We see such features in the profile close to the star. If the temperature and pressure 
are high enough (T > 1340 K, p > 10 -5 bar) so that even refractory solids would sublimate, the estimation 
of the elemental abundance breaks down. Inside this radius r in , listed in table [2] no condensates will form. 
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Figure 1: From top to bottom: Surface density profiles, temperature profiles and pressure profiles in the midplane, the m5 
disk in the left column, the mlO disk in the right column at transition time. Each plot shows the two-dimensional model (solid 
line), the Chambers model (dashed line) and the self-similar solution model (dotted line). The adopted model parameters are 
given in table [2] 
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Figure 2: The weight percentage chemical abundance of some of the most frequent elements in solid species as a function of 
distance from the star. The left panel gives the two-dimensional model (2D) in the case m5, and the right panel gives the 
Chambers disk model (Ch) in the case mlO. The elements arc from top to bottom at 4 AU: O, Fe, Si, S and Al. For example: 
in the 2D case, a solid particle condensing at 1 AU is composed by weight of 30 % O, 30 % Fe, 20 % Si and 20 % is comprised 
of all other elements. 



The self-similar solution model and the two-dimensional model in the mlO disk fail to provide condensates 
across the whole annulus a — [0.5 AU, 4 AU] where planetesimals in the dynamical simulations arc initially 
located. 

In general, close to the star, only refractory species condense out of the solar nebula. The abundance of 
more volatile elements increases with distance from the star. Hence, a hot disk model generates a disk 
with a relatively high abundance of refractory elements. On the other hand, cooler disks provide higher 
abundances of volatile elements. 



3.3. Source regions in the N-body simulations 

In the dynamical simulations, the orbits of the planetesimals and protoplanets change due to various mech- 
anisms, such as close encounters or secular resonances. Gas drag damps the eccentricities of the smallest 
planetesimals and moves them inwards. If they do not grow fast enough and do not merge with other bodies, 
they do not decouple from the gas. Hence, the planetesimals at small initial semi-major axis fall into the 
star. On the other hand, larger bodies start to migrate inward due to type I migration as long as the gas 
disk has a significant surface density. The gas dissipation time scale is of the order of few Myr. As a result, 
most of the material initially located inside 1 AU is lost in almost every simulation. 

The region of initial locations of the planetesimals that end up embodied in a final planet is what we call 
the planet's source region. In order to study the source regions of individual planets, we group planets first 
according to their semi-major axis. We divide the planets in four groups (a, b, c and d, ordered by increasing 
semi-major axis), each containing the same total mass (and thus the same number of initial planetesimals 
(see figure [3|). To study the source region as a function of the planet's final mass in each of these groups, 
we divide each group a, b, c and d in two sub-groups containing the same number of initial planetesimals: 
a group of high mass planets and a group of low mass planets. With the statistics of 4 orbital groups and 
4x2 mass sub-groups, the dependence of the source region on the planet's semi-major axis and final mass 
can be studied. A study of the dependence on the simulation initial conditions would be desirable but the 
number of simulations is too small. A planet's source region is characterised by the cumulative distribution 
function of planetesimals as a function of radius from which it is built. We compare the source regions of 
any two groups using the two-sided Kolmogorov-Smirnov-test. The estimated p-value gives the probability 
that the two distributions of the initial planetesimal locations are drawn from the same parent distribution. 
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Figure 3: Left: The initial distribution of all planetesimals that end up in one of the planets in four orbital groups indicated 
by different colors (a: black, b: red, c: blue, d: green). The lower panel shows the distribution of all planets in radial distance, 
their diameter represents their masses. Right: The distribution of the planetesimals of the small mass (cyan) and high mass 
(magenta) sub-groups in the four orbital groups a-d. The lower panel in each histogram gives the distribution of the orbital 
group that is considered. 

If p is small, it is unlikely that the source regions coincide. 

First, we explore how the source regions differ depending on the final orbital radii of the planets. For 
each planet in a given orbital group, we gather the initial positions of its building blocks. Figure [3] gives 
the histograms of the locations of all the initial planetesimals that end up within the planets of the four 
different orbital groups. Planets close to the star have a wide and flat source region (group a), they form 
from material that was also closer to the star initially. Planets at larger distances (groups b, c and d) all 
have steeper and narrower source regions. For these groups, the source regions arise entirely from beyond 
a radial distance of 1 AU. This is mainly because almost all material initially within this location ends up 
in the star. For example, the black histogram (group a) indicates that more than 200 planetesimals (« 5% 
of the total mass of the group) initially located within 1.0 AU end up in the planets in group a, whereas 
according to the green histogram, not more than 10 planetesimals end up in planets of group d. KS-tests 
reveal that the groups do not share the same initial distribution as all p- values are below 10~ 6 . This means 
that we should not expect that planets at different semi-major axes have the same source region. The inner 
most group a stands out most significantly from the others (see figure [3]) . 

Next, we study how the source regions differ depending on the final mass of the planets. The planets in each 
orbital group are split into two sub groups according to their mass. The KS-test shows that the sub-group 
containing low mass planets and the sub-group containing high mass planets do not have the same parent 
distribution in all of the orbital groups. A trend is visible: low mass planets have a source region that is 
steeper and narrower, whereas massive planets have a wider and flatter source region (see figure (31). 

Often, the planets are located much closer to the star than the average radii of the initial planetesimals 
from which they were built. This is due to orbital migration. Most of the bulk composition of the planets 
is contributed by the solids that condense out of the solar nebula at 1 AU< r < 4AU, except some of the 
very close-in planets. 

3.4- Combining disk chemistry and N-body simulations 

In the next step, the calculations of the chemical abundance in the disk model and the dynamical simulations 
are combined. According to our transition criterion, the disk model dependent transition time tdisk estimated 
in section |3.1| is the time in the disk model evolution at which the chemical composition of solids in the 
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1.58 


S, % 


0.24 


1.62 


0.89 


2.2 


Ca, % 


1.18 


1.61 


1.6 


1.33 


C, ppm 


5.1 


468 


44 


2960 


N, ppm 


0.046 


4.3 


0.59 


180 


Na, ppm 


200 


1390 


2450 


5770 


P, ppm 


390 


1860 


1200 


1100 


Ti, ppm 


630 


850 


800 


650 


Cr ppm 


7180 


4060 


3400 


3680 



Table 3: Planetary a bundances in wt.% or weight ppm. References: ^Morgan and Andersj ( |i980fr , ^Kargel and Lewis| (1993 
^Lodders and Fegley|(l997) . 



disk represents the composition of the planetesimals at the beginning of the dynamical simulations. The 
initial radial distribution of the planetesimals that end up in a given planet gives the final bulk chemical 
composition of this planet. The element abundance profiles in the disk and the diversity in the source 
regions suggest that different planets have different bulk chemical compositions. All planetesimals from the 
dynamical simulations that contribute to the mass of any planet originate from beyond the condensation 
line ri n at ttrans- There the composition of all planets can be fully determined. The element abundance in 
the planets for each disk model is shown in the online supplementary data. 



The simulations by Morishima et al. (2010) are designed to reproduce the inner planets of the Solar System, 



thus, we now compare the element bulk abundance of the simulated planets with the Solar System terrestrial 
planets. 



A normalized abundance Nx of element X relative to Si can be obtained from (34) via 



N x,t P = 



( mL%ufx \ 
\wt.%oiSi J 



sp 



/ wt. %ofx \ 
{wt.%ofSi) tp 



(39) 



where the wt.% is based on simulated planets (subscript sp) or literature values of terrestrial planets (sub- 



script tp) in our Solar System. The bulk chemical compositions were taken from Morgan and Anders| ( 1980 ) , 
Kargel and Lewis (1993), and Lodders and Fegley (1997), see table [3] One should keep in mind that these 
values are partly based themselves on simulations and that uncertainties lead to absolute errors up to the 
order of 25% on these abundances (Bond et al. 2010a). These uncertainties are discussed in the caveats 
section 14.11 



Next, we search for a simulation outcome that best reproduces the composition of the inner Solar System 
planets. We choose a planetary system including four planets out of the mlO simulations and a similar 
system out of the m5 simulations. Since all planets have to form in the same simulation and the number of 
simulations is limited, we focus on systems that match the constraints of our Solar System best according to 



Morishima et al. (2010), keeping in mind that differences in mass and semi-major axis have an effect on the 



source region. EJS 2-1-5 (short hand for r gas = 2Myr, p = 1, m5 (see section 2.3)) and EJS 3-1-10 are two 
systems that fulfill many of the constraints, although the masses of the planets are very poorly reproduced: 
Venus and Earth are too small and the Mars and Mercury analogs are too massive. The source region of 
close in planets is sensitive to the planets mass and the large Mercury mass in the EJS 2-1-5 simulation 
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Mercury 


Venus 


Earth 


Mars 


EJS 2-1-5 


mass [m$] 


0.52 


0.45 


0.79 


0.51 




a [AU] 


0.31 


0.59 


0.94 


1.39 


EJS 3-1-10 


mass [m®] 


0.09 


0.23 


0.34 


0.27 




a [AU] 


0.67 


0.89 


1.26 


1.54 


Solar System 


mass [mj,] 


0.06 


0.81 


1.00 


0.11 




a [AU] 


0.38 


0.72 


1.00 


1.52 



Table 4: Parameters of simulated planetary systems that represent the Solar System, a is the semi-major axis of the planet. 
The Solar System values are taken from |Kargel and Lewis] jl993| . 



might affect the result. The large Mars mass is a general problem that is often observed in dynamical 
simulations that reproduce the inner Solar System, recently discussed in |Walsh et al. (2011). The masses 



and semi major axis of the chosen simulated planets are shown in table |4j 

In figure [4j the normalized abundance iVx,tp * s snown f° r each element and all Solar System planets. Both 
systems fail in reproducing in detail the bulk chemical composition of the Solar System. The most refractory 
elements (Al, Ti and Ca) are underestimated in Mercury, Venus and Earth analogs. There are two reasons 
for this: Gas drag and type-I migration leads to the loss of refractory rich planetesimals that form in the inner 
part of the disk. Hence, they can not be embedded in the planets. On the other hand, if the temperature of 
the disk is high enough so that refractory solids dominate the solid abundance not only in the very center 
of the disk, more planetesimals are dominated by refractory material. The Chambers model in the mlO 
case predicts the highest overall temperature and the normalized abundances of refractory material for the 
Earth and Mars analogs are significantly increased. The most volatile element abundances (Na and S) are 
overestimated. This discrepancy is higher if the disk is too cold. In this case volatile rich planetesimals 
dominate the disk and the normalized bulk abundance becomes too high in the case of all four planets. 
Except in the case of S, the composition of Mars is reproduced by most of the disk models. On the other 
hand, Mercury is not very well reproduced in any of the models. 

Note that the differences in bulk composition of the planets we are considering is small. In both simulations, 
the abundances differ by up to 50 % in the case of refractory elements if the Chambers-model is used. All 
the other element abundances nearly coincide. The fact that the variation in bulk chemical composition in 
the simulated planets is in general smaller than in the Solar System rocky planets is also demonstrated in 
figure [5| 

Terrestrial planet formation and especially the chaotic growth phase when many giant impacts take place is 
a stochastic process. This might result in the formation of unique planets that reproduce the bulk chemical 
composition of certain Solar System planets very well, but to reproduce the bulk composition of all Solar 
System planets in a single simulation might be a very rare situation. Hence, we study the distribution 
of the element abundance of all planets at once to get an overview as to how extreme the composition of 
planets in those simulations can be. As source regions of different planets can differ significantly, we expect 
variations in the element abundances depending on the underlying disk model. In figure [5] the distribution 
of the elemental bulk chemical compositions of planets is shown for each element in the case of the most 
extreme disk models normalized to the Earth values. One model is the two-dimensional model in the case 
m5. This is the model that gives the lowest temperature across the disk and a lot of species condense 
across the planetesimal annulus. Thus, it results the most inhomogencous clement abundance profile. The 
other model is the Chambers model in the mlO case where the resulting element abundance profiles are the 
most featureless since this disk has the highest temperature. Although these models are applied one two 
different samples of planets, comparing their element abundances reveals the most extreme variations that 
are possible in our approach. 

In the two-dimensional model, the relative abundances of the simulated planets are mostly located in a 
narrow range. In the case of refractory elements, the abundances are a little lower than the Earth value and 
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Figure 4: The normalized abundance of elements for Mercury, Venus, Earth and Mars analogues, from bottom to top. The 
elements are listed in order of increasing volatility from left to right. Two-dimensional model (solid line), Chambers model 
(dashed line) and self similar solution model (dotted line) disks are taken into account. Left column: EJS 2-1-5, right column: 
EJS 3-1-10. 
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they agree with Mars. There is a variation in O, especially for small semi-major axes. In the case of the 
volatile elements Na and S, the bulk chemical composition is very high relative to the Earth values, up to a 
factor of 9 in the case of sulfur. Some planets with small semi-major axis tend to have a smaller abundance 
of these volatile elements. 

In the Chambers model, there is a variation in the relative abundances, up to a factor of 5 in the case of 
sulfur. Concerning refractory elements, some planets still group along a narrow range at the same position 
as in the two-dimensional model. They form a lower bound of the distribution. Some planets have much 
higher clement bulk abundances and agree with the relative abundances of Mercury. Beside some outliers, 
these planets have small semi-major axes. The volatile element abundances are reproduced in a similar 
manner than before, although the variation towards lower relative abundances is a little larger. 

The bulk chemical compositions are much more sensitive to the different source regions than in the cooler 
two-dimensional model. This can be explained using the example of sulfur: in the two-dimensional model 
FeS condenses around 1 AU, hence most planets incorporate planetesimals with a similar weight percentage 
of sulfur. The few outliers result from the difference in the source regions. In the Chambers model, condensed 
sulfur is only present outside 2.5 AU. Hence, depending on the source region, a larger variation in the bulk 
composition is possible. 

Another example is given by the refractory element Al. In the hot disk, Al is the dominant element inside 
the condensation radius of more volatile solids at around 1.3 AU (see figure [2]). Hence, planetesimals that 
are locate initially inside this radius have a wt.% that is higher than any planetesimals that form in the 
colder disk. This results in the higher variation. 

In both models, the differences between the simulated planets and the Solar System planets are significant, 
especially in the volatile regime and in the case of Mercury. 

Geochemical ratios allow us to quantify the depletion in volatile elements (e.g. Na/Si) with respect to 
the enrichment in refractory elements (e.g. Al/Mg). The ratios in the Solar photosphere (table [I]) are 
Na/Si = 0.045 and Al/Si = 0.072. For the Earth, according to table [3| Na/Si = 0.017 and Al/Si = 0.098. 
This reflect the volatile depletion and refractory enrichment in the Earth relative to Solar System bulk 
compositon (Palme and O'Neill 2003 Rubie et al. 2011[ ). Figure [5j which gives the ratio normalized by the 



reference Earth ratio, shows that the depletion in volatilcs is not achieved in any of the simulated planets. 
Most planets have a Solar photosphere ratio. On the other hand, the suitable enrichment in refractories is 
provided by some planets in the warm disk model. 



3.5. Water delivery 



In the dynamical simulations by O'Brien et al. (2006), the EJS and EEJS simulations result in a water content 



of Earth-analogs that is very low, since with high giant planet eccentricities the supply of planetesimals from 
the icy asteroid belt region are suppressed. With the same giant planet orbits, in the simulations created by 



Morishima et al. (2010), a larger fraction of the planets originates from the outer region of the planetesimal 



annulus, which should result in a overabundance of water. 

Since water is one of the most important molecules in the formation and evolution of life, we focus briefly on 
the H abundance in the simulated planets to quantify the above statement by Morishima et al. (2010). We 



keep in mind that the condensation calculations will not be valid at temperatures where hydrated species 
form. Nevertheless, the abundance of H within a planet in our model will provide a benchmark on the 
abundance of material from the outer edge of the planetesimal disk. 

The H abundance is mainly controlled by the condensation of serpentine. All disks are too hot for water ice 
to condense. Serpentine condenses out of the solar nebula in the outer region of the planetesimal annulus 
within 4 AU before itrans depending on the disk model. This edge is located inside 4 AU only in the self- 
similar solution model and in the two-dimensional model. In figure [6j the number of planets comprising a 
wt.% of H is shown. In general, different disk models lead to a completely different wt.% of H since its inner 
edge of condensation is close to the outer edge of the dynamical simulations. 
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Figure 5: The element abundance of all planets. The ordering of the elements increase in volatility from top to bottom. The 
horizontal axis gives the abundance of element X of a planet relative to the abundance of Si, see ( |39| l normalized by Earth 
values. The vertical axis gives the semi-major axis of the planet. The diameter of the black dots represents the mass of the 
simulated planets. The red markers indicate location of Mercury (left triangle), Venus (square) and Mars (right triangle), the 
dashed line highlights the Earth abundance. The left column represents m5 with a two-dimensional disk model (2D), this is 
the model that predicts the lowest disk temperature. The right column represents mlO with a Chambers disk model (Ch), this 
is the model that predicts the highest disk temperature. 



19 



2D m5 



1.0 0.1 0.2 0.3 0.4 0.5 
Wt% of H 




1.00 0.01 0.02 0.03 0.04 0.05 
wt% of H 



Figure 6: The wt.% of H of all planets for the two-dimensional (2D) and the self similar solution (Se) disk models for m5 disks. 
All the other models give no condensation of H rich material in the disk. 



3.6. Extreme disk chemistry 



In Bond et al. (2010b), the bulk chemical composition of hypothetical terrestrial planets in extrasolar 



planetary systems is simulated. The element abundance of stars can be significantly different compared to 
the abundance in our solar nebula and giant planets in other systems can have completely different masses 
and orbits than Jupiter and Saturn. These circumstances could lead to the formation of rather extreme 
objects such as carbon-rich planets. 



We explore one extreme case by adopting element abundances of HD 4203 used in Bond et al. (2010b). A 



recent article by Fortney (2012) pointed out that the derived stellar C/O ratios could be overestimated. 
Nevertheless, we use the same disk models and dynamical simulations as before. If we focus on the most 
extreme disk models again, figure [7] reveals that the wt.% of C becomes extremely high for small radial 
distances. The main forms of carbon at these distances are C, TiC and SiC. Considering that most of the 
planetesimals are initially located inside ~ 1 AU, planetesimals that will end up in a planet do not have 
concentrations higher than 40 wt.%. From this, we can expect that a wt.% of 40 might be the highest bulk 
abundance that can be achived in the planets. The C-abundance of the planets normalized to the Earth 
value is shown in figure [8] Note that around one quarter of all planets do not contain C and are not plotted. 
This also holds in the case of the two-dimensional disk model. In both disk models, there is a wide spread 
in the C-abundance which can be up to a factor 10 4 times higher than the Earth C-abundance. Hence, the 



wt.% of C in the planets range from a few to 40 percent. This verifies results by Bond et al. (2010b I. 



4. Caveats, Implications and Future Work 

In this section, we focus on the major caveats and implications that result from our approach concerning 
the three main goals we introduced at the beginning: First, the estimation and interpretation of abundance 
profiles in the circumstellar disk. Second, the interpretation of the source regions of the simulated planets 
in the dynamical simulations. Finally, the combination of abundance profiles and dynamical simulations, 
the dependence of the resulting bulk chemical compositions on the initial assumptions and an application 
to the Solar System. 



4-1. Caveats 

Disk abundance profiles. Strong assumptions and simplifications are required when estimating chemical 
abundance profiles. There are a lot of uncertainties in describing the disk physics correctly, for example the 
source of viscosity in the disk and the opacity. Solid migration is not taken into account in the disk models, 
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Figure 7: These figures show the weight percentage chemical abundance of the most frequent elements in solid species as a 
function of radius in the case of extreme initial abundances. The models are from left to right: Chambers with mlO and 
two-dimensional with m5. The thick line shows a very high abundance of C close to the star. At 4 AU, the elements are from 
top to bottom: O, Fe, Si, Mg, C. 
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Figure 8: The C-abundance of all planets in the case of a extreme solar nebula. Details are explained in the caption of figure [5] 
(no reference values of Mercury, Venus and Mars are shown here). Note the different scaling of the abscissa since the variations 
are much stronger. The left panel represents mlO with a Chambers disk model (Ch). The right panel represents m5 with a 
two-dimensional disk model (2D). 
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although it is expected to occur in the actual circumstellar disks (Garaud 20071. Here, although we used 



very simple models, they give rise to a broad range of temperature, pressure and surface-density profiles. 
The self-similar solution model we constructed is not fully self-consistent in its derivation. Since we did not 
use it in the further study, this caveat is not important. 

All calculations performed with the HSC Chemistry code are made under the assumption of chemical 
equilibrium. Each solid species is treated as a separate phase, and each is considered a pure substance. 
Therefore, only the most stable species appear as condensates. Actually, the formation of solids can be 
controlled by disequilibrium condensation, at least partially. In such a case, once a solid has formed, it does 
not interact anymore with the system. This may occur if condensation takes place fast enough and could have 
a significant effect on the estimated bulk composition ( |Bond et al. 2010a]). Other disequilibrium effects arc 



evaporation or photo-dissociation of gaseous molecules (Ebel 20061. At low temperatures, the assumption 



of chemical equilibrium is not valid anymore due to lack of diffusion in the condensed solids. Hence, we can 
not make reliable predictions for the volatile content in solids and the condensation of hydrated silicates like 
serpentine, the only solid containing H in our calculations, might be not justified. However, condensation 
of some hydrated material can not be excluded at the very outer edge of the planetesimal belt. Beside this, 
the radial abundance of biologically important elements like H, C and N can not be estimated correctly with 
our approach. C and N rich species do not condense in any of the models. They would condense out outside 
of 4 AU, were the disk is cold enough. In addition, the formation of water out of H and O depends on its 
physical environment and the presence of H on a planet need not indicate the presence of H 2 0. 

The transition from the disk models to the dynamical simulations contains uncertainties. Our transition 
approach is based on matching the initial surface density of solids at 1 AU by adapting the time of the 
disk model evolution. The slope of the dynamical simulation surface density is exactly reproduced in the 
self-similar model. The other models provide a flatter profile and the surface density profile deviates from 



the initial conditions in Morishima et al. (2010 1. We did not study the dependence of the abundance profiles 



on the transition criterion systematically. If the surface density at 1 AU is changed, as it is the case for 
the two different solid disk masses in the dynamical simulations, the temperature and pressure profiles and 
thus the element abundance in the disk change significantly. A similar effect occurs when the position of 



the normalization is changed. If the transition time is left as a free parameter (as in Bond et al. (2010a)) 



the abundance profile and thus the bulk compositions will change in a similar manner as if the planetesimal 
disk mass is changed. If the transition takes place later, the temperature is lower which corresponds to a 
smaller disk mass. 

Another caveat is the gas-to-dust ratio we adopted to estimate the surface density of solids in our disk 
models. The Morishima et al. (2010) simulations use a gas-to-dust ratio which is up to 3 times higher. 



This high ratio would result in a very high surface density at ttrans- The temperature and pressure in the 
disks would be very high, which moves r- m to a larger radial distance. Thus, no solids would be present in 
large parts of the circumstellar disk at itrans- This would be inconsistent with the initial condition of the 
dynamical simulations. Therefore, we claimed that a gas-to-dust ratio similar to the ISM is more reasonable. 
Although this results in an inconsistency with the initial conditions of the dynamical simulations, a smaller 
gas density in the dynamical simulations would only affect the gas-drag and type I migration. The strength 



of the effect of this type of migration is still a matter of debate (Schlaufman et al. 2009 Paardekooper and 



Papaloizou 2009) 



Both caveats concerning the transition mentioned above could be eliminated by adding a more sophisticated 
disk model to the dynamical simulations which could directly provide temperature and pressure profiles and 
such a transition would be eliminated. 



Source regions of planets. Observations reveal that there is a very large variety in the configuration of 
extrasolar systems. Due to limitations in computational resources, the simulations by Morishima et al. ( 2010 1 
were constructed in order to reproduce the constraints given by the Solar System and their outcomes do not 
represent general planetary systems. The presence of Jupiter and Saturn truncates the initial planetesimal 
disk at 4 AU and no icy bodies beyond this limit are taken into account. Schlaufman et al. (2009) have 
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shown that the type-I migration rate obtained from linear theory decreases by an order of magnitude, if 
for example the ice line is taken into account. These recent results are not included in the dynamical 
simulations. In addition, the analytic expression for type-I migration used by Morishima et al. (20101 



does not account for the sharp surface density transition which may occur at the inner edge of the gas 
disk, where the disk is truncated by e.g. magnetosphcric accretion. Since these strong density gradients 



can act as protoplanet traps (Masset 2006 Hasegawa and Pudritz 2011), they can prevent early forming 



protoplanets from falling into the star. Moreover, at the pressure maximum at the inner gas disk edge, gas 
drag will accumulate planctcsimals and the edge can also act as a planctcsimal trap. These effects might 
have significant consequences on the source regions and the final location of the planets. 



Since only one simulation per set of initial conditions was performed by Morishima et al. (20101, we are not 



able to do a broad study of the dependence of the source regions on the initial conditions of the dynamical 
simulations. The late stage of formation is a stochastic process and the source regions of two planets in one 
group might naturally differ significantly between each other. A much larger set of simulations is necessary 
to explore stochasticity and parameter space. 



Combining chemistry and dynamics. Our transition approach is based on matching the initial surface den- 
sity of solids at 1 AU by adapting the time of the disk model evolution. This should not imply that all 
the planetesimals in the dynamical simulations formed instantaneously at their initial positions. Yet, all 
planetesimals are allocated the composition of solids in the disk only at the time of transition. Actually, the 
formation and growth of planetesimals out of solids might significantly depend on the local density, temper- 
ature and pressure and does occur inhomogeneously along radial direction in the disk. A more sophisticated 
chemical model could provide condensation rates of solids, so that one could keep track of the solid fraction 
at every radius and integrate them over time to improve our snapshot-like approach. 

The inner edge r; n of the disk of solids that condenses at itrans lies inside the inner edge of the initial 
planetesimal annulus in the mlO simulations for all disk models. Thus, no solids can be allocated to 
these planetesimals. However, in all of the dynamical simulations, all the planetesimals initially located 
inside r- m are lost to the star due to gas drag and type I migration and the bulk composition of the final 
planets can be completely estimated in every case. Nevertheless, two caveats remain: a smaller migration 



rate (Paardckooper and Papaloizou 2009 Schlaufman et al. 2009) could prevent some of the innermost 



planetesimals from migrating into the central star. In addition, the gravity of the planctcsimals initially 
located inside r ln gravitational interact with their counterparts on larger radial distance. Hence, they can 
not simply be ignored. 



According to Albarede et al. ( 2009 1 , the inner part of the Solar System was very hot as long as gas was 



present and any solids that condensed were depleted in volatiles relative to the Solar photosphere. Our disk 
models do not predict such high temperatures at transition time and the expected volatile depletion in the 
case of e.g. the Earth is not achieved. However, at the inner edge of the planetesimal annulus, highly volatile 
depleted planetesimals form but most of them are lost into the star. As previously mentioned, modifications 
in the migration mechanism might result in a higher fraction of volatile depleted planetesimals that end 



up within planets. Following |Albarede et al. (2009J, the local planetesimals were dry and the water has to 
come from asteroids and comets from the Jupiter-Saturn region or comets from the trans-Uranian region 



(Morbidelli et al. 2000). However, water- vapor absorption on silicate grains can be a mechanism to form 



water-rich planetesimals in the hot regions of the disk ( Muralidharan et al. 2008) 



The reference values for the bulk composition of the Solar System rocky planets were not reproduced in 
detail. Since these values are themselves based to some degree on simulations and condensation calculations, 
they do not necessarily provide real constraints. The bulk composition of Earth (Kargel and Lewis 1993) 



is based on Best Bulk Silicate Earth (BBSE) abundance values, which are estimated using several methods 
that have their own uncertainties. Knowing the BBSE and the volatility trends for the elements, the core 
abundances and the bulk chemical composition of the Earth was estimated. Uncertainties in the BBSE 
values, the volatility trends and the ratio of core mass to total mass result in relative errors from 5-10% 



for most elements and up to 25% for volatiles. In Lodders and Fegley (1997), element abundances for Mars 
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are based on meteroite samples and an oxygen isotope model. The uncertainties result in a relative error of 
10%. In Morgan and Anders (1980), since the data is very limited, the compositions of Venus and Mercury 



are determined fully theoretically. In their method, Solar gas cools under the assumption of equilibrium 
condensation. The composition of a body is determined by the amount of early condesates which were most 

The uncertainties arising from this method are 



1978) 



likely embodied in the final planet (Morgan et al. 
not quantified. In summary, the fact that we can reproduce parts of the characteristica of the chemical 
compositions of the planets does not imply definitly that our model is reproducing the Solar System planets. 
The reference values should be seen as a guideline. 

There are some further caveats resulting from the dynamical simulations: the limited number of simulations 
does not provide any planetary system that accurately mimics the inner Solar System. Since terrestrial 
planet formation is a stochastic process, a larger set of simulations could provide better candidates. The 
simulations adopt perfect sticking if two bodies merge, and no decrease of volatiles is taken into account. 



More sophisticated accretion models ( Stewart and Leinhardt 2009 1 propose a collision outcome depending 
on the collision parameters concerning energy loss and fragmentation. It ranges from a perfect accretion 
event to complete disruption of impactor or target. The energy release of these events is clearly high enough 
that volatile elements can be lost. This can be a possible explanation of the overabundance of volatile 



elements (Bond et al. 2010a), although there are other additional processes that can result in the depletion 
in volatiles in the solar nebula (Davis 20061. Finally, we note that the rather extreme composition of 



Mercury, e.g. its massive core, can be explained by a collision that striped off most of Mercury's mantle and 
crust. No such event can take place in the Morishima et al. (2010) simulations. 



4-2. Implications 

Disk abundance profiles. Different initial conditions in the dynamical simulations result in different abun- 
dance profiles in the disk. Since the temperature and pressure in a disk is directly related to the surface 
density, the choice of the initial disk mass and the profile of the planetesimals in the dynamical simulations 
controls the distributions of condensed solids in the disk. Relatively cold disks shelter a smaller variety of 
solids and especially smaller abundances of volatiles like H. A massive disk of planetesimals such as the mlO 
simulations implies a hot disk where only a few solids will condense close to the center of the disk and no 
H is embodied in solids inside 4 AU. A less massive disk like the m5 disk can result in the condensation of 
the entire range of considered elements across the disk. The initial disk mass has an equally large or even 
larger effect on the resulting abundance profile than the choice of the disk model. 



Source region of planets. Comparing the source regions of the planets in the different groups reveals that 
most of them do not coincide. The distributions arc sometimes more peaked or more flat, but the width 
of the source region is always very similar. Due to migration and secular resonances, the final planets are 
located in the inner part of their source region, some of them are even found outside of their source region. 
The innermost planets tend to have source regions that are flatter. A trend is visible in the planet mass 
dependence of source regions: Massive planets have a flat and wide source region while small planets have 
a steep source region which becomes steeper at larger radial distance. This implies that massive planets 
consist of elements that condensed from regions that spanned a larger range of physical conditions than 
their less massive counterparts. 



Combining chemistry and dynamics. The bulk chemical composition of planets produced in the dynamical 
simulations do not reproduce the bulk chemical composition of the rocky Solar System planets in detail. 
We studied two realizations of a four-planet system: In both cases, the effect of different disk models or the 
initial composition was not dramatic and the two dynamical simulations also gave similar results. Thus, it 
may be difficult to reproduce the diversity in the reference bulk chemical compositions of our Solar System 
in a single simulation. 



In Bond et al. (2010a), the transition time is a free parameter and a different parameter space is covered 



than in our study. We believe that our transition criterion is more self-consistent, but both methods result 
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in a broad range of planetary compositions. Similar to Bond et al. (2010a), we overestimated the abundance 



of volatile elements, implying that volatile loss can not be ignored. Keeping in mind the differences in the 



dynamical simulations by Morishima et al. (20101 and O'Brien et al. (20061, we found that it is the choice 



of the disk model (a cold disk) that controls the amount of wet planetesimals that are accreted (rather than 
the orbits of the Jovian planets) . 

A comparison of the bulk chemical composition of all simulated planets reveals the effect of the choice of the 
disk model. If the initial disk mass of the dynamical simulations is small (m5) and the coldest disk model is 
chosen (two dimensional model) , most of the simulated planets share the same bulk composition. Variations 
are only visible in O, Na, and S for planets with small semi-major axis. In this scenario, Solar System values 
are almost never achieved. This means that the dynamics are not significantly affecting the abundances of 
this model. On the other hand, a massive disk (mlO) with a hot disk model (Chambers model) provides 
large variations in the most refractory and volatile elements, since dynamics are important in this scenario. 
Again, small semi-major axis planets tend to vary more than others. 

A similar trend is visible in the case of planet mass, especially for S, since more massive planets vary more 
in bulk composition than smaller planets. A difference was also observed for the source regions of large 
mass and small mass planets. Massive planets tend to be comprised of a more uniform distribution of 
planetesimals that formed all across the disk. Therefore they tend also to vary more in bulk composition 
than small mass planets, which are comprised of planetesimals from the outer disk region. A larger set 
of dynamical simulations is needed to quantify this difference. Note that such correlations between mass 
and bulk composition are not visible in the extreme abundance plot for C (figure [8]) . The C abundance 
dominates over a larger region than for example Al in the case of a Solar nebula composition and thus, the 
different source regions of small and large mass planets do not play an important role. 

Both the self-similar solution model and the two-dimensional disk model in case of the less massive disk, give 
rise to the condensation of H rich material inside the planetesimal region. Water is an important component 
in the formation and evolution of life as we know it. If predictions about water on habitable planets in such 
simulations are made, the choice of the disk model and the initial disk mass have a significant effect on the 
amount of H that is delivered to the planets. However, condensation calculations fail in this temperature 
regime and perhaps the planetesimals form dry and water has to be delivered from outside the planetesimal 
belt. 



4-3. Future work 

Given all of the uncertainties and assumptions, the main caveat in this work is the limited number of 
dynamical simulations that were available. A new code currently under development at the University of 
Zurich that will be able to perform terrestrial planet formation simulations using CPUs and a much larger 
number of simulations is planned to explore stochasticity and a broader parameter space. With a more 
sophisticated disk model in the dynamical simulations, estimating the bulk composition of planetesimals 
should be possible more accurately and self-consistently. A large set of dynamical simulations will also 
provide a good starting point for statistical studies on the source region of planets. Future studies should also 
include a more detailed treatment of the condensation sequence which would provide a better understanding 
of the evolution of solid abundances in the disk. 



5. Conclusions 

The combination of chemical disk models and dynamical simulations opens a new avenue for exploring planet 
formation models and to interpret the history of the terrestrial planets ( |Bond et al. 2010a). The output of 
the coupled disk model, its chemistry and the dynamical simulations leads directly to estimates of the bulk 
composition of the simulated planets. We explored the sensitivity of the radial element abundance trends 
in the resulting planets to the disk model assumptions, the initial conditions of the dynamical simulations 
and initial composition of the solar nebula. Our main conclusions are: 
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Bond et al. (2010b I concluded that their estimated bulk chemical composition are in excellent agree- 



ment with observed planetary values, indicating that the models in the dynamical simulations of 



O'Brien et al. (20061 work properly in reproducing the bulk composition of the inner Solar System. 
In our study, these compositions are not reproduced in detail, resulting from the different dynamical 
simulations and the more restrictive but self-consistent transition from the gas disk to the planetes- 
imals. The largest discrepancies are in the abundances of volatiles and refractories and in the bulk 
composition of a Mercury analogue. 

• The source regions of planets are unique. Large mass planets form from a broader region of the initial 
disk whilst lower mass planets tend to comprise of planetesimals which are initially located in the 
outer region of the disk. 

• The element abundance profiles change significantly if different disk models or if different initial disk 
masses of planetesimals are used. Some elements do not condense at all if the disk model predicts 
a high temperature along the disk. This dependence has a significant effect on the estimated bulk 
chemical compositions of the planets. 

• In a cold disk model, the variety in the bulk compositions is small and mixing is not important. In a 
hot disk model, the dynamics lead to a large variety in the bulk compositions of planets in the case of 
the most refractory and most volatile elements. The variation tends to be larger for massive planets 
and for those with small semi-major axis. 
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Online supplementary material 

Tables [5] to |l0| show the bulk chemical composition of all planets in the simulations for different disk models 
and disk masses. They are published as online material. 
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